/*

File name		: figureE1.do

Description		: Grid search of the parameter γ (ITT version in Appendix)
									
Output File		: figureE1.png
				
*/

set more off, perm
capture restore

*** 001 Grid Search of parameter gamma ***
use "data/blood_donation/blood_donation.dta", clear

	// Merge rainfall data
	merge m:1 entnahmedatum using "data/rainfall/rainfallZH.dta"
	*, keep(match)

	// Rename kalendertag daily_rainfall
	rename daily_rainfall_ZH daily_rainfall
	gen any_rain = 1 if daily_rainfall > 0 & daily_rainfall != .
	replace any_rain = 0 if daily_rainfall == 0
	sort spnr entnahmedatum
	gen any_rain_plus1 = any_rain[_n+1] if spnr == spnr[_n+1]

	sort spnr period
	sum don if period == 1
	sum don if period == 2 & don[_n-1] == 1 & spnr == spnr[_n-1]
	sum don if period == 2 & don[_n-1] == 0 & spnr == spnr[_n-1]
	sum don if period == 3 & don[_n-1] == 1 & spnr == spnr[_n-1]
	sum don if period == 3 & don[_n-1] == 0 & spnr == spnr[_n-1]
	sum don if period == 4 & don[_n-1] == 1 & spnr == spnr[_n-1]
	sum don if period == 4 & don[_n-1] == 0 & spnr == spnr[_n-1]

	/* Gen new variable that is 1=donated in period 1, 2=did not donate in period 1, 
	 3 = donated in period2, 4=did not donate in period 2 */
	sort spnr period
	gen treatment_bydonation1=1 if period == 1 & don == 1
	replace treatment_bydonation1=0 if period == 1 & don == 0
	replace treatment_bydonation1=treatment_bydonation1[_n-1] if spnr == spnr[_n-1]
	gen treatment_bydonation2=1 if period == 2 & don == 1
	replace treatment_bydonation2=0 if period == 2 & don == 0
	replace treatment_bydonation2=treatment_bydonation2[_n-1] ///
	if spnr == spnr[_n-1] & treatment_bydonation2[_n-1]!=.
	replace treatment_bydonation2=treatment_bydonation2[_n+1] ///
	if spnr == spnr[_n+1] & treatment_bydonation2[_n+1]!=.

// Define grid search program
global within_type .
capture program drop within2
program within2

	local idfe_yes $idfe_yes

	capture drop loglik1 loglik2
	capture drop loglikk1 loglikk2
	
	matrix ll = (.,.)
	matrix lll = (.,.)
	
	forvalues i= 0.1(0.01)0.9 {
		
		capture drop phonecall_vector
		gen phonecall_vector = phonecall + `i' * phonecall_minus1 + ///
		`i'^2 * phonecall_minus2 + `i'^3 * phonecall_minus3
			
			if `idfe_yes' == 0 {
				qui reghdfe don phonecall_vector ///
				male age bltype_dummy2 bltype_dummy3, ///
				absorb(c_num)
			}
			else{
				qui reghdfe don phonecall_vector, ///
				absorb(c_num spnr)
			}
						
			matrix ll[1,1] = `i'
			matrix ll[1,2] = `e(rss)'
		
		if `i' == 0.1{
			matrix loglik = ll
		} 
		else{
		 	matrix loglik = (loglik \ ll)  
		 } 
	}
	
	*max
	svmat loglik
	sort loglik2 loglik1
	sum loglik1 if _n==1&loglik1[1]>0.1&loglik1[1]<0.9
	
	local max=`r(mean)'
	local lb=`max'-0.01
	local ub=`max'+0.01
	
	forvalues j= `lb'(0.001)`ub' {
		
		capture drop phonecall_vector
		gen phonecall_vector = phonecall+`j' * phonecall_minus1 + ///
		`j'^2*phonecall_minus2 + `j'^3*phonecall_minus3
			
			if `idfe_yes' == 0 {
				qui reghdfe don phonecall_vector ///
				male age bltype_dummy2 bltype_dummy3, ///
				absorb(c_num)
			}
			else{
				qui reghdfe don phonecall_vector, ///
				absorb(c_num spnr)
			}
			
			matrix lll[1,1]=`j'
			matrix lll[1,2]=`e(rss)'
			
		if `j'==`lb'{
			matrix loglikk = lll
		} 
		else{
			matrix loglikk = (loglikk \ lll) 
		}  
	}
	
	*max
	svmat loglikk
	sort loglikk2 loglikk1
	matrix max = .
	sum loglikk1 if _n == 1&loglikk1[1]>0.1&loglikk1[1]<0.9
	local rounded = round(`r(mean)', .001)
	capture matrix max = `rounded'

end

// Run within2 first without individual fixed effects
global idfe_yes 0
	within2

	gen loglik2_fig=loglik2
	gen loglik1_fig=loglik1

// Then run within2 first with individual fixed effects
global idfe_yes 1
	within2

	graph twoway (scatter loglik2_fig loglik1_fig, color(navy) yaxis(2) ///
	ylabel(758.5 (0.5)760.5, axis(2)) yscale(axis(2) range(758.5 760.5) ///
	titlegap(3)) text(427.7 0.1 "Blood drive FE" "& individual FE", ///
	place(n) color(navy)) text(427.37 0.1 "Blood drive FE", place(n) color(navy)) ///
	text(427.77 0.47 "{&gamma}=0.470""  (0.129)", place(n) color(maroon)) ///
	text(427.89 0.47 "{&beta}{subscript:1}=0.093""  (0.016)", place(n) color(black))) ///
	(scatter loglik2 loglik1, yaxis(1) color(navy) ///
	graphregion(color(white)) xscale(titlegap(3)) ylabel(426.7 (0.3) 427.9, axis(1)) ///
	yscale(axis(1) range(426.7 427.9) titlegap(3)) bgcolor(white) ///
	ytitle("Residual sum of squares", axis(1)) ytitle("Residual sum of squares", axis(2)) ///
	xtitle("Habit formation parameter") legend(off) ///
	text(427.77 0.655 "{&gamma}=0.655""  (0.191)", place(n) color(maroon)) ///
	text(427.89 0.655 "{&beta}{subscript:1}=0.125""  (0.023)", place(n) color(black))) ///
	(scatteri 426.7 0.655 427.75 0.655, color(maroon) yaxis(1) recast(line) lpattern(dash)) ///
	(scatteri 758.5 0.470 760.25 0.470, color(maroon) yaxis(2) recast(line) lpattern(dash))

	graph export output/figures/figureE1.png, replace width(2048)	
